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We  present  the  results  of  an  investigation  of  the  Prony— Lanczos  (P-L) 
method  {14,38]  and  the  power  method  [39]  for  simple  computation  of 
approximations  to  a  few  eigenvectors  and  eigenvalues  of  a  Eermitian  matrix. 

We  are  motivated  by  realization  of  high-resolution  signal  processing  in  an 
integrated  circuit.  The  computational  speeds  of  the  above  methods  are 
analyzed.  They  are  completely  dependent  on  the  speed  of  a  matrix-vector 
product  operation.  If  only  a  few  eigenvalues  or  eigenvectors  are  needed,  the 
suggested  methods  can  substitute  for  the  slower  methods  of  the  LINPACK  or 
EISPACK  subroutine  libraries.  The  accuracies  of  the  suggested  methods  are 
evaluated  using  matrices  formed  from  simulated  data  consisting  of  two 
sinusoids  plus  gaussian  noise.  Comparisons  are  made  with  the  corresponding 
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I.  Introduction 


Ve  are  motivated  by  the  use  of  eigenvector  decompositions  of  data 
matrices  or  estimated  covariance  matrices  for  detection  of  signals  in  noise 
and  for  estimation  of  signal  parameters.  This  has  evolved  from  early  work  of 
Liggett  [1]  and  Owsley  [2],  to  adaptive-array-detection  improvements  of 
Tufts  and  Kirsteins  [3,33]  and  high-resolution  parameter  estimators  of 
Cantoni  and  Godara  [4],  Bienvenu  and  Kopp  [5],  Owsley  [6],  Schmidt  [21]  and 
Tufts  and  Eumaresan  [7,32]. 

Principal  component  analysis,  using  principal  eigenvalues  and 
eigenvectors  of  a  matrix,  was  initiated  by  Karl  Pearson  (1901)  [8],  and 
Frisch  (1929)  [9]  in  the  problem  of  fitting  a  line,  a  plane  or  in  general  a 
subspace  to  a  scatter  of  points  in  a  higher  dimensional  space.  Eckart  and 
Young  [34]  presented  the  use  of  singular  value  decomposition  for  finding 
low-rank  approximations  to  rectangular  matrices.  C.R.  Rao  examined  the 
applications  of  principal  component  analysis  [10].  Eigenvector  analysis  is 
also  used  in  image  processing  to  provide  efficient  representations  of 
pictures  [11].  Recently,  principal  component  analysis  ras  been  coupled  with 
the  Wigner  mixed  time-frequency  signal  representation  to  perform  a  variety 
of  signal  processing  operations  [28,30,31]. 

Linear  Prediction  techniques  for  estimation  of  signal  par ameters , which 
are  modern  variants  of  Prony's  method,  can  be  improved  using  eigenvector 
decomposition  [7].  Prony's  method  is  a  simple  procedure  for  determining  the 
values  of  parameters  of  a  linear  combination  of  exponential  functions.  Now 
'Fiuny's  method*  is  usually  ttken  to  mean  the  least  squares  extension  of  the 
method  as  presented  by  Hildebrand  [13].  The  errors  in  signal  parameters 
which  are  estimated  by  Prony's  method  can  be  very  large  [14].  If  the  data 
is  composed  of  undamped  sinusoids,  the  forward  and  backward  prediction 
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equations  and  a  prediction  order  larger  than  the  number  of  signal  components 
can  be  used  simultaneously  as  idvocated  by  Nutall  [22],  Clrych  and  Clayton 
[23],  and  Lang  and  McClellan  [24].  Tufts  and  Kumaresan  have  shown  how  one 
can  improve  such  methods  of  parameter  estimation  by  going  through  a 
preprocessing  step  before  application  of  Prony's  method  [7,15,16,17]. 
The  measured  data  matrix  or  the  matrix  of  estimated  covariances  is  replaced 
by  a  matrix  of  rank  M.  which  is  the  best  least  squares  approximation  to  the 
given  matrix.  If  there  is  no  prior  information  about  the  value  of  M,  it  is 
estimated  from  the  data  using  singular  value  decompositon  (SVD) . 

The  eigenvalue  problem  [37]  is  one  area  where  extensive  research  has 
been  done  and  well  established  algorithms  are  available  in  highly  optimized 
mathematical  libraries  such  as  LINPACK  and  EISPACK  [40]  .The  computational 
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complexity  of  these  algorithms  is  of  order  0(N  )  where  N  is  the  size  of  the 
matrix.  They  solve  for  the  complete  set  of  eigenvalues  and  eigenvectors 
of  the  matrix  even  if  the  problem  requires  only  a  small  subset  of  them  to  be 
computed.  For  the  above  appl icat ions , only  a  few  principal  eigenvectors  and 
eigenvalues  are  needed.  Hence, we  would  like  to  use  a  method  which  uses  this 
specialization  to  reduce  the  computations. 

Tufts  and  Kumaresan  [29,32,33]  have  suggested  procedures  for 
improving  Prony's  method  without  computation  of  eigenvectors.  These  appear 
to  perform  about  the  same  as  the  more  complicated  approaches  which  use 
eigenvalue  and  eigenvector  decomposition.  The  approach  in  [29]  is  based  on 
the  results  of  Hocking  and  Leslie  for  efficient  selection  of  a  best  subset 
[25].  The  approach  of  [32]  and  [33]  is  based  on  the  simple  computations 
which  result  from  using  the  longest  possible  prediction  interval. 


Fere  we  investigate  two  different  approaches  to  achieving  S\T>-like 
^nprovenent  to  Frony's  method  without  the  coopu t a t ; cna  1  cost  of  actually 
computing  the  SM)  or  computing  all  eigenvectors  and  eigenvalues.  The  idea 
is  to  calculate  the  f ew , ne ce s sary  eigenvalues  and  eigenvectors  using  the 
power  method  [39]  and  a  method  of  Lanczos  [14].  Our  derivation  of  Lanczos' 
method  stresses  the  connection  with  Prony's  method  .  The  methods  are 
analyzed  and  their  amounts  of  computation  are  calculated.  Simulations  are 
performed  and  results  are  compared  to  the  singular  value  decomposition 
method  in  LINPACX. 


I I .  The  Prony-Lanc zo s  Method 

Let  us  assume  that  we  start  with  a  given  square .Hermit ian  matrix  A  for 
which  we  want  to  compute  the  principal  eigenvectors  and  eigenvalues.  For 
examples,  this  could  be  either  the  true  underlying,  population  covariance 
matrix  or  the  estimated  covariance  matrix  [36]  from  spatial  or  temporal 
data.  Let  us  also  define  the  eigenvectors  and  eigenvalues  associated  with 
the  matrix  A  (dimension  A=n). 

Aui  =  ,  i  =  1,2, ...,n  (1) 

where  =  0  ,  i  /  j 

u^*'uj  =  1  ,  i  =  j  ,  that  is  u^  are  orthonormal  vectors.  (2) 

The  asterisk  is  used  to  denote  a  complex  conjugate  transpose. 

The  characteristic  polynomial  associated  with  the  matrix  A  is  given  by 
det(A  -  AI)  =  0  (3) 


Expanding  the  determinant  we  have  the  polynomial  equation 


Pn-1 


n- 1 


+  p  =0 
^  n 


(4) 


and  the  roots  of  this  polynomial  will  give  us  the  eigenvalues  of  the 

matrix.  We  briefly  summarize  the  procedure  for  obtaining  the  eigenvalues  Xj 


based  on  the  Lanczos  'power  suns"  as  presented  in  [14].  We  shall  show  that 
the  eigenvalues  can  then  be  obtained  from  the  power  sums  by  Frony's  method 


[13]. 

Let  us  select  a  starting  vector  b^.  We  assume  that  the  starting  vector 
^  has  a  non-zero  projection  on  the  eigenvectors  of  the  matrix  A 
corresponding  to  the  eigenvalues  that  we  want  to  compute. 

We  then  analyze  the  vector  in  the  reference  system  of  the  vectors 
{u ,  which  are  the  set  of  orthonormal  eigenvectors  of  the  matrix  A: 

-o  ~  '^1  5i  ^^2  H2  "^n  5n 

where 


T  •  =  u  ■  b 
1  -1  1-0 


(6) 


Hence,  using  equation  (1) , successive  vectors  formed  by  premultiplications 
of  by  powers  of  the  matrix  A  can  be  represented  as  follows  : 


-1  “  ^^o  ~  ■'l  ^1  ll  “'2  ^2  22  ^  "^n  ^n  2n 


b2  =  =  Abj 


tj  uj  +  T2  22  ^ 


.  T_  X_  U_ 
n  n  — n 


(7) 


^k+1  =  ^o  =  Ab^  =  Ti  Xj^  +  ^2  >^2’'  22  +  '  •  •  ^  ^n  ^n’'  2n 
Let  us  form  the  set  of  basic  scalars: 

‘^i+k  ^  ^i  ^k  ""  ^k  ^i 
Then  we  shall  have: 

‘^k  =  hlP  ^1^  ^  h2l^  X2^  ^  +  1tJ2  x^k  =  (9) 

which  were  called  by  Lanczos  the  "weighted  power  sums"  [14]  . 

The  problem  of  obtaining  X^'s  from  the  Cj^'s  is  the  "problem  of  weighted 
moments*  [14]  .  That  is  the  problem  of  Prony  [12]  and  the  old  and  modern 
versions  of  Prony's  method  can  be  used  to  estimate  the  X^'s. 

The  prediction-error-filter  equations  of  Prony's  method  can  be  written 
as  follows: 
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eg  +  C,g,  +...+C  "  =0 

o°o  11  a-1  u-1  n 

Cig  +  c-g  +...+cg  ,+c^,  =0 

1  o  21  n  n-1  n+1 


or  in  matrix  form. 


(10a) 


eg  +  c  -g.  +  .  .  .  +  c-  ,g  ,  +  c-  =  0 
n  o  a+1  1  2n”l  n**!  2n 


c  ■  s  =  0 


(10b) 


A  non-zero  solution  is  possible  if  the  determinant  of  C  is  zero. 

From  the  theory  of  Prony's  method  [13] 

g(X.)  =  +  g^_l  +  .  .  .  +  gl  Xi  +  =  0  (11) 

hence  the  polynomial  coefficient  vector  g  is  also  orthogonal  to  the  vector 
2  k  X 

(1  Xi  Xi  ...  where  X  ^'s  are  the  eigenvalues  of  the  matrix  A. 

Lanezos  noticed  that  Prony's  method  can  be  simplified  if  we  i^ubstitute 
the  sequence  {1  X^  Xj^  ...  X^”}  for  a  row  of  the  matrix  C  to  form  a  matrix 
C'.  If  we  replace  the  matrix  C  by  C'  in  (10b),  the  non-zero  vector  ^  is 
still  a  solution,  because  of  (11).  Hence  the  determinant  of  C'  must  be  zero. 


det  C’ 


X.^  ,  .  .  X.° 

1  1 


...  c 

2  n 


=  p'(X.)  =  0 


c  c  . 

n-1  n 


Hence,  the  X;'s  can  be  obtained  directly  by  finding  the  zeros  of  the 


polynomial  p'(z).  That  is,  Lanezos  showed  that  it  is  not  necessary  to  first 


solve  equations  (10)  for  the  prediction-error- f ilter  coefficients. 

Thus,  in  the  absence  of  noise,  we  know  that  entering  the  weighted  pov/er 
sums  c^  of  (8)  in  equation  (12)  and  finding  the  roots  of  the  resulting 

polynomial  will  provide  us  with  accurate  estimates  of  the  true  eigenvalues 
of  the  covariance  matrix  A.  Note  also,  that  equation  (12)  can  be 
reduced  to  a  2*^*^  order  equation  involving  only  c^,  c-^,  cj  and  still 

provide  us  with  accurate  solutions  for  our  problem  of  estimating  one  or  two 
s inuso ids . 

Now,  if  our  data  is  composed  of  one  or  two  complex  sinusoids,  then  the 
(LxL)  covariance  matrix  elements  will  be  also  one  sinusoid  or  a  sum  of  two 
sinusoids,  hence  the  rank  of  the  matrix  will  be  one  or  two, respectively. 
The  eigen-decomposition  of  the  matrix  will  show  that  it  has  only  one  or  two 
non-zero  eigenvalues  and  hence  it  can  be  characterized  by  a  linear 

combination  of  one  or  two  eigenvectors,  corresponding  to  the  principal  non¬ 
zero  eigenvalues.  In  Appendix  A  it  is  shown  that  these  eigenvectors  can  be 
expressed  as  a  linear  combination  of  complex  sinusoids  which  have 
frequencies  equal  to  these  of  the  sinusoids  composing  the  data. 

Now, suppose  that  we  have  accurately  determined  a  few  e  igenvalues , say 
two,Xj  and  (nxn)  matrix  A.  We  wish  to  determine  the 

corresponding  eigenvectors.  Two  concepts  are  used  :  (a)  premultiplication  of 
a  vector  by  the  matrix  (  A-X^I  )  removes  the  i'"^  eigenvector  component  of 
that  vector  and  (b)  if  a  vector  ,  to  a  good  approximation, consists  only  of  M 
eigenvector  components  , then  removing  (M-1)  of  these  components  leaves 
one, isolated  eigenvector  component. 

Let  us  consider  the  special  case  of  a  rank  two  matrix  : 

^  ~  ^1-lSl*  ^252^2*  (13) 
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From  eouations 


'  5  )  and  (13 )  we  have : 


Abo  =  ‘'•2^252 

Then, our  pre 1 iminary .unnorna 1 ized  estimates  of  the  two  principal 
eigenvectors  are  : 


Ui'  =  (A-X2l)Abo  =  (A-X2I )  (■C2^llil'‘"'-2^2-2^  “ 

_  2  ,  .  2  2 

-  "^2X2  22  +  '^2^2  -2  ~  ^1^1^251  ~  "^2^2  -2  ~ 

~  ^ iX 1 ( Xi~X2 ) Uf 


(15) 


And  similarly  for  the  second  eigenvector  estimate  we  have 


—2*  ~  2^2 ^ ^2"”^  1  ^ ^2  (16) 


Normalizing  the  eigenvectors  u^'  (i=l,2)  we  can  write  (15)  and  (16)  as 
22'  =  e'^^l  2i  '  angle  of  r2X2(X2"X2)  (17) 

22'  =  52  ;  ^2“  angle  of  T2''-2 ^^2'*^1  ^  (18) 


In  general  .given  the  required  eigenvalues  from  the  earlier  Prony 
calcnlat ion.we  estimate  an  unnormalized  eigenvector  from  the  formula  : 


uj^'  =  T  (  A-Xjl  )  Abo 
i^k 


(19) 


where  the  number  of  factors  in  ;  he  product  depends  on  the  number  of 
significant  eigenvector  components  in  Ab^. 

Finally,  a  few  comments  should  be  made  on  the  selection  of  the  starting 
vector  bg.  Our  sole  assumption  until  now  has  been  that  b^  has  a  non-zero 

projection  on  some  eigenvector  of  A  that  we  want  to  compute.  A  good  b^ 

vector  would  have  to  be  biased  in  favor  Oi  e  principal  eigenvectors.  t7e 
have  found  that  the  Fourier  vector  provides  a  very  good  selection  for  b^. 
This  vector  will  have  its  fundamental  frequency  computed  from  the  maximum 
peak  of  the  DFT  data  spectrum.  Very  frequently  in  signal  processing 
applications  the  data  is  preprocessed  through  a  DFT  step  for  a  coarse 
analysis.  This  is  a  valuable  ’■onus  for  our  method  to  use  the  available 

information  for  further  processing. 

Ill .  The  Power  Method 

Suppose  A  is  a  Hermit ian  (nxn)  matrix.  The  SVD  theorem  [37]  states 

that  A  can  be  written  as: 

A  =  D  •  S  •  (20) 

where  U  is  a  unitary  matrix  and  S  is  a  matrix  consisting  of  real  only 
diagonal  elements  [37]. 

The  power  method  computes  the  dominating  singular  vectors  one  at  a  time 
and  is  based  on  solving  the  equation: 

su  =  An  (21) 
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for  the  singular  vector  n  and  the  singular  value  s. 
an  iterative  scheme  to  solve  (21).  (Ve  instead  sugge 
using  an  appropriate  starting  vector  b^: 


The 


t  a 


power  method  uses 
two-step  solution 


Bi  =  A  /  ||A  b^ll 


(22) 


The  singular  value  is  chosen  to  be: 


M  =  llA 


(23) 


In  order  to  obtain  the  next  singular  vector,  the  estimated  singular  plane 
X 

(ulUi  )  is  removed  from  A  using  the  following  deflation  procedure  [37]; 


A' 


=  A  -  sj  Si  u^^ 


(24) 


and  the  procedure  is  repeated  with  matrix  A  to  yield  S2»U2‘ 

The  selection  of  b^  is  very  important  and  the  Fourier  vector  provides  a 
very  good  estimate.  This  preprocessing  step  can  be  implemented  in  VLSI  very 
efficiently  using  summation-by-parts  [28]  or  the  Fast  Hartley  Transform 
[42,43]  methods.  A  necessary  thing  r-^quired  to  implement  the  power  method 
is  a  circuit  capable  of  computing  matrix  vector  products  of  the  form  An. 
But  the  rounding  errors  associated  with  it  are  always  worrisome  limiting  the 
usefulness  of  the  power  method.  For  this  reason  we  propose  to  use  the 
permuted  difference  coefficients  (PDC)  algorithm  [26,27]  coupled  with  the 
known  Fourier  vector  to  perform  the  above  operation  with  high  accuracy  and 
no  round-off  errors.  A  VLSI  implementation  for  the  PDC  algorithm  can  be 
easily  realized  using  a  random  access  memory  (RAM)  toghether  with  a  read- 
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only-memory  (ROM)  where  the  original  Fourier  coefficients  and  the  subsequent 
reordered  coefficients  addresses  are  stored. 


IV.  Operation  coxint 

In  this  section  we  calculate  the  total  operations  needed  for  the 
singular  value  decomposition  (LINPACK),  the  Prony— Lanczos  method  and  the 
Power  method. 

(1).  The  matrix  eigenvalue  problem  has  been  solved  in  both  LINPACK  and 
EISPACK  mathematical  libraries.  The  LINPACK  SVD  routine  is  presented  here. 
The  solution  can  be  divided  in  three  steps:  reduction  to  bidiagonal 
form,  initialization  of  the  right  and  left  unitary  matrices  U  and  V  and  the 
iterative  reduction  to  diagonal  form. 

The  reduction  to  bidiagonal  form  has  the  following  floating  point 
multiplication  count  (for  a  square  NxN  matrix): 

2[  N^  -  N^/3] 

Approximately  the  same  number  of  additions  are  required. 

In  the  second  step  the  amount  “^f  work  involved  when  only  the  right-hand 
side  matrix  V  is  computed, is: 

2N^/3 

floating  point  multiplies  and  approximately  the  same  number  of  additions. 

In  the  last  step  rotations  are  used  to  reduce  the  bidiagonal  matrix  to 
diagonal  form.  Thus  the  amount  of  work  depends  on  the  total  number  of 
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rotations  needed.  If  this  number  is  r,  then  we  have  the  following 
multiplication  counts: 


4Nr 


The  number  r  is  quite  difficult  to  estimate.  There  exists  an  upper  bound  for 
r , 


r  <  sN^/2 


where  s  is  the  maximum  number  of  iterations  required  to  reduce  a 
superdiagonal  element  as  to  be  considered  zero  by  the  convergence  criterion. 
Hence  the  total  operation  count  for  the  LINPACK  SVD  solution  is: 

2N^+4Nr  <  2N^(s+l)  flops 

where  by  the  term  "flop*  we  denote  a  floating  point  multiply-add  operation. 

(2).  The  Prony-Lanczos  method  is  entirely  dependent  on  the  speed  of  a 
matrix-vector  product  operation  .  For  a  rank  two  square  matrix  of  size  N  we 
shall  have: 

The  matrix-vector  multiplications  to  determine  the  vectors  b^^  involve 
N  floating  point  multiplications  and  (N-1)  floating  point  additions  per  row 
for  a  total  of  : 


flops 


2 

(  2N  flops  for  the  two  vectors  bj.b2  )•  The  scalar  weights  c^  .i=0,l,2,3 
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require  vector-vector  inner  products  for  a  count  of  N  multiplications  and 


(N-1)  additions  per  vreight  .  Therefore  the  total  is; 

4N  flops 

The  computation  of  the  eigenvalues  from  the  (second  order)  determinant 
condition  involves  12  flops  and  one  square  root  calculation.  Finally, the 
eigenvector  computation  requires  N  flops  for  each  vector  for  a  total  of  2N 
flops . 

Hence  the  total  operation  count  for  the  Prony—Lanczos  procedure 
requires : 


(2N^+6N+12)  flops  +  1  square  root 

The  above  computations  do  not  include  the  work  required  to  select  the 
starting  vector  using  a  DFT  analysis.  In  this  case, assuming  a  data 

sequence  zero  padded  to  M  points, we  shall  have: 

Mlcg2M  flops 

plus  (M-1)  additions  for  the  determination  of  the  maximum  spectral  peak. 

(3).  The  power  method  computes  the  dominating  eigenvalues  and 
eigenvectors  one  pair  at  a  time  .  The  second  pair  will  be  computed  following 
a  deflation  of  A.  In  general,  the  number  of  iteration  steps  depend  on  the 
convergence  criterion  severity  .  We  instead  claim  that  two-steps  are 
generally  enoug'  o  provide  sufficient  accuracy.  The  Fourier  vector  is  again 
selected  as  the  starting  vector  b^ . 
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Pne  first  e  igenva  lue /e  igenvec  tor  pair  requires  2N  +2N'  flops.  The 

2  2 

deflation  step  requires  N  flops  and  N  floating  point  additions. 

Hence  (for  a  rank  two  matrix)  the  power  method  requires  a  total  of 

5N^+4N  flops 

2 

plus  N  floating  point  addtions. 

V.  Simulation  results 

Let  us  assume  that  we  have  a  data  sequence  which  is  composed  of 

uniformly  spaced  samples  of  two  closely  spaced  complex  sinusoids  in  white 
noise.  We  shall  follow  the  methods  described  earlier  in  section  II  i  III  to 
calculate  the  principal  eigenvalues  and  eigenvectors. 

The  data  sequence  is  given  by  the  equation 

x(n)  =  exp(j2nfj^n  +  9)^)  +  exp(j2nf2n  + 

with  =  0.52Hz,  f2  =  0.5Hz  and  for  n=l,2,...,25 

Here,  25  data  samples  are  used  and  the  phase  difference  is  A0  =  n/2 
computed  at  the  middle  of  the  data  set,  effectively  reducing  the  signal-to- 
noise  ratio  in  that  region,  thereby  representing  the  worst  case  that  can  be 
encountered.  The  frequency  separation  is  less  than  the  reciprocal  of  the 

observation  time.  The  data  is  zero  padded  to  m=128  points  and  then  the 
maximum  peak  of  the  DFT  is  computed  to  yield  the  frequency  of  the  Fourier 

vector.  This  vector  will  be  used  as  a  starting  eigenvector  for  the  P-L  and 

Power  methods  later. 

We  construct  the  forward  plus  backward  augmented  covariance  matrix  A  of 
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size  (21x21)  .  Its  effective  rank  is  two.  The  SVD  routine  .the  P-L  method 
and  the  Power  method  are  employed  to  solve  for  the  eigenvalues  and 
eigenvectors  (eigenpairs)  of  the  matrix.  The  P-L  method  and  the  Power 
method  compute  only  the  two  principal  eigenpairs.  The  mean  values  and 
standard  deviations  of  the  eigenvalue  estimates  are  given  in  Table  I  for  an 
ensemble  of  500  experiments.  The  performance  of  the  P-L  and  Power  methods  is 
almost  identical  to  the  SVD  (LINPACK)  method  for  the  first  eigenvalue 
estimates.  At  high  SNR  the  second  eigenvalue  mean  and  standard  deviation 
estimate  obtained  from  the  P-L  method  is  biased  with  respect  to  the 
noiseless  SVD  results.  However  .at  low  SNR  the  eigenvalue  statistics  ire 
closer  to  the  noiseless  SVD  results  than  the  other  two  methods. 

Table  II  presents  the  statistics  of  the  distances  of  the  P-L  and  Power 
methods  eigenvectors  from  those  of  the  SVD  method.  The  distance  is  the 
inverse  cosine  of  the  angle  between  the  subspaces  spanned  by  the  estimated 
eigenvectors  [41].  The  results  show  that  for  the  first  eigenvector  the  P-L 
estimate  of  the  mean  is  less  biased  (about  one  order  of  magnitude)  than  the 
Power  method,  whereas  for  the  second  eigenvector  estimates  they  perform  the 
same.  This  shows  that  these  vectors  span  virtually  the  same  subspace  as  the 
vectors  computed  from  the  SVD  method.  The  eigenvector  estimates  were  also 
compared  to  the  signal  eigenvectors  and  the  distances  were  computed  as 
above.  The  results  show  that  at  high  SNR  the  eigenvector  spanned  subspaces 
have  a  greater  distance  from  the  signal  subspace  than  the  SVD  subspace.  At 
low  SNR  the  distance  is  reduced  and  the  second  eigenvector  statistics  are 
closer  to  the  signal  eigenvector  than  the  S\'D  eigenvector. 

Table  III  shows  the  CPU  time  reiiuired  to  compute  the 
eigenvalues/eigenvectors  pairs  for  these  methods.  The  P-L  method  is  faster 
than  the  SVD  by  the  order  of  the  size  of  the  covariance  matrix,  which  here 
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is  21.  This  roughly  agrees  with  the  theoretical  operation  count  we  presented 
in  section  IV.  It  is  almost  twice  as  fast  as  the  Power  method.  Inclusion  of 


the  FFT  computation  in  these  two  methods  will  offset  some  of  their  speed 
advantage  over  the  SVD  .  Nevertheless  .the  P-L  method  is  again  about  one 
order  of  magnitude  faster  than  the  SVD  method  and  the  Power  method  a  little 
more  than  half  that  (6  times  faster). 

The  frequencies  f^  are  then  obtained  from  the  eigenvectors  of  the 
estimated  covariance  matrix  by  the  T-K  method  [7].  For  both  estimates  of 
the  mean  and  standard  deviation  ,as  presented  in  Table  IV, all  three  methods 
perform  similarly  down  to  15  db.  At  Odb  the  P-L  method  yields  slightly 
better  statistics  than  the  other  two  methods. 


VI .  Cone  fusion 

Two  methods, the  Prony^Lanezos  method  and  the  Power  method  are  proposed 
for  simple  computation  of  approximations  to  a  few  eigenvectors  and 
eigenvalues  of  a  Hermitian  matrix.  The  computational  speeds  of  these  methods 
were  analyzed.  The  accuracies  of  the  proposed  methods  were  evaluated  using 
covariance  matrices  from  data  consisting  of  two  sinusoids  in  a  gaussian 
noise  environment.  Comparisons  were  made  with  the  corresponding  eigenvectors 
and  eigenvalues  obtained  using  the  LINPACK  mathematical  library.  The 
suggested  methods  can  substitute  for  the  slower  method  of  LIKPACK  if  a  few 
eigenvalues  or  eigenvectors  are  needed. 
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In  this  appendix  we  derive  the  eigenvalues  and  eigenvectors  of  the 
covariance  matrix  R  for  the  case  of  one  and  two  sinusoids. 

One  Complex  Sinusoid  Case: 

The  data  sequence  is  modelled  by: 

jUj^n 

y(n)  =  a^e  ,  n  =  1,2,. ...N 

The  covariance  values  of  y(n)  are: 


N 


1  ^  • 

r^(i,j)  =  2  y  (n-i)y(n-j)  i,  j=l,2,  . .  .  ,L 

n=L+l 


(A.l) 


Writing  the  covariance  matrix  R  explicitly  in  terms  of  the  signal,  we  have; 


|a^|^  "•'‘^1 


,2 


l»ll  « 


i^r 


Ki 


2  -ju)j(L-2) 


j(Oj(L-l) 


|aj|  e 


Kr 


(A.2) 


We  can  diagonalize  R  by  an  orthogonal  matrix  U  resulting  in  the  following 


equa  t ion ; 


The  eigenvalnes  of  R  which  occur  along  the  diagonal  elementsnof  the  above 
equat ion, sat isfy  the  following  equation; 


L 

^  X,  =  tr(R)  =  L  ja^]^  (A. 4) 

i=l 


But  the  covariance  matrix  R  is  of  rank=l,  since  it  has  only  one  linearly 
independent  row  (or  column).  The  rest  are  obtained  by  multiplying  by  a 

+  •  i 

- 

constant  number  (e  ). 

Then  the  eigenvector  corresponding  to  the  eigenvalue  Xj  =  L  is; 


ju)  2j(a  j(L-l)(i) 
Uj=c^(le  e  e  ) 

ir 


can 


since  it  annihilates  every  row  of  the  matrix  (R-X.^1).  The  constant  c^ 
be  determined  from  the  fact  that  the  matrix  D  is  orthonormal,  hence: 


which  yields; 


1 


Hence  finally; 


e 


j (L-l)u 

e  ) 


and  this  is  a  Fonrier  vector  with  fundamental  frequency 
Two  Complex  Sinusoids  Case: 

The  data  sequences  is  modelled  by; 


(A. 5) 


IS 


'(n)  =  a. 


J‘a2n 


n  =  1.2, 


Tlie  covariance  estimates  are  given  by  the  expression: 


r  (k.m) 

yy 


N 

^  y  {n-k)y(n-m) 
n=L+l 


I  ® 


+  \a.^\  e  +  +  a^a^  ,k.m= 

where : 


V 


1 


1 

N-L 


1 


n=L+l 


j  (n-m)- jw^  (n-k) 

e 


N 


n=L+l 


juj  (n-m)- (n-k) 
e 


Rewriting  the  matrix  R,  we  have; 


R  = 


M 


2 


(Lx2) (2xL) 


where  : 


(A.  6) 


^  (A. 7) 


(A. 8) 
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Ml  = 

M2  = 

M^i  1^  ii 
til* 

^  £2 

e 

22  -* 

1^2!^  22  ^  ^*  2l^ 

2l  = 

[1  e>i 

_  2  j<t) 

e  J  1 

,j(L-l)u,^jT 

22  = 

[1  e->^2 

,j(L-l)a>^jT 

and 


^1*2  V  -j 

*  =-N=L -  ®  2  ® 

n=L+l 


If  Uj  is  an  eigenvector  of  R  corresponding  to  eigenvalue  Xj,  then: 


^1^221  “  ^l2l  (A. 9) 

Premultiplying  by  H21  have: 

■^2^1^221  “  ^2  2l  (A.  10) 

Thus  Xj  is  also  the  eigenvalue  of  Vl2^1  corresponding  eigenvector  is: 

II  =  *^221  (A.  11) 

Premultiplying  (A. 10)  again  by  Mj 

^1^2^12i  ~  ^l^l2l  (A. 12) 

and  comparing  (A. 10)  with  (A. 12) 

2l  =  %Il  (A.  13) 


Thus  we  can  find  the  eigenvalues  and  eigenvectors  of  R  by  working  with  the 
matrix  M2M2  which  is  of  order  2.  Hence: 
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The  eigenvalaes  Xj  and  X2  are  found  to  be: 

Xj  =  l/2CL|ai|^  +  L|a2p  +  2Re{xg)  +  ((L|ai|^  +  L|a2p+  (A. 16) 

2Re(xg})2-4(L2  -  |g |^ ( [a^ p  1*2 p- |x ^ 

X2  =  l/2{L|aiP  Lh2p  2Re{xg]  -  ((L|aj^P  +  L|a2p  +  2Re{xg}^)  - 
4(l2  -  |8p)(|.i|2|.jp  - 

whe  re 


Re{xg}  =  Re 


*  N 

{a^a2  ]  cos(— 

2(N-L) 


-2  )A«i) 


A(i)  s  in 


(A. 17) 


Note  that  a  column  of  the  adjoint  of  (M2Mj-XjI)  gives  the  eigenvector  Vj  of 

M2M1. 
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Adj 


(LjaJ^  -  ig)  - 


* 

-  Lx 


+  X  g  )- 


T 

Therefore  the  eigenvector  v^  is  :  Zl  ~  ^^21^ 

or  vj  =  [(L|aj|^  +  xg)-Xj^  -  |aj^  |^g-Li3^ 

Now  the  eigenvector  Uj  of  R  corresponding  to  Xj  is: 

Zl  "  ”lZl 


and  hence. 


Zl  =  ’^11  0^1  P  Zl  ^Z2^  '^2lUa2p  Z2  **  Zl^' 
a  linear  combination  of  the  Fourier  vectors  £2  and  £2* 

Similarly,  the  eigenvector  £2  of  ®  corresponding  to  X2 

Z2  =  "^ll'M^lPzi  ^^21' ^  l®2pZ2'^**Zl^ 

T 

where  v^ '  =  ''2l'^ 

and 


is : 


vii'  =  (L|a2|^  +  r*g*)-X2 

'^21'  =  -  |»2P  «  " 

The  rest  of  the  eigenvalues  of  R  are  zero  and  the  corresponding 
are  not  unique. 


(A. 18) 

1 


(A. 19) 


(A. 20) 


(A. 21) 


(A. 22) 

e igenvectors 
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• 

SNTl 

SVD 

P-L 

PM 

mean=  22.0357 

22.0126 

22.0174 

st.dcv=  0 

0 

0 

22.0636 

22.0423 

22.0353 

30 

0.2652 

0.2655 

0.2642 

22.0341 

22.3182 

22.2957 

15 

1.4927 

1.4936 

1.4892 

28.8561 

28.5285 

28.5425 

0 

8.6489 

8.7576 

8.7477 

Eigenvalue  estimate 


I 


1.7107 

0.5741 

1.7131 

03 

0 

0 

0 

1.7162 

0.7504 

1.7199 

30 

0.0497 

0.3777 

0.0498 

1.8634 

1.0327 

1,8677 

15 

0.2856 

0.5357 

0.2884 

10,5379 

1.6797 

7.3859 

0 

2.7981 

1.2625 

2.6301 

L 

Eigenvalue  estimate  X2 


■ 


TABLE 


SNR 

S’V'D.F-L 

SVD,PM 

mean  =  0 

0.4770  -4 

OD 

st.dev  =  0 

0 

0.698C  -5 

0.5819  -4 

30 

0.1938  -4 

0.7575  -4 

0.2775  -4 

0.1991  -3 

15 

0.5746  -4 

0.1606  -3 

0.5305  -2 

0.5932  -2 

0 

0.1019  -1 

0.1035  -1 

1- 

First  Eigenvector  Distances 


C9 

0.3917  -4 

0 

0.6169  -4 

0 

0,1744  -3 

0.3283  -3 

30 

0.1162  -3 

0.2850  -3 

0.4618  -2 

0.2055  -2 

15 

0.3682  -2 

0.1110  -2 

0.8243  -1 

0.7146  -1 

0 

0.2614  -1 

0.2938  -1 

- 

Second  Eigenvector  Distances 


TABLE  II 
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SNR 

S\T) 

FFT 

P-L 

PM 

-35 

0.30472 

+  5 

0.14050 

+4 

0.15835 

+4 

0.27610 

-r4 

30 

0.24391 

+5 

0.14119 

+4 

0.15859 

+4 

0.27793 

+4 

15 

0.24538 

+5 

0.14065 

+4 

0,15877 

+4 

0.27506 

+4 

0 

0.25819 

+5 

0.14029 

+4 

0.15874 

+4 

0.27568 

+4 

Computational  Cost 

(  measured  in  time  units  t^, where  1  tg=  26,04166  p.sec  ) 


TABLE 


III 


SNR 

SVD 

P-L 

PM 

mean=  0.5000 

0.5000 

0.5000 

CD 

st.dev=  0 

0 

0 

0.4999 

0.4999 

0.4999 

30 

0.0013 

0.0013 

0.0013 

0.4961 

0.4952 

0.4962 

15 

0.0157 

0.0137 

0.0154 

0.4331 

0.4620 

0.4551 

0 

0.1334 

0.0898 

0.1082 

. 

Frequency  Estimate  fj 


CO 

0.5200 

0 

0.5200 

0 

0.5200 

0 

30 

0.5201 

0.0013 

0.5201 

0.0013 

0.5201 

0.0013 

0.5251 

0.5249 

0.5251 

15 

0.0190 

0.0141 

0.0190 

0.5717 

0.5613 

0.5642 

0 

0.1184 

0.0893 

0.0980 

. 

Frequency  Estimate  £2 


TABLE  TV 
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